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Abstract 

In this work, some of the NLO QCD corrections for pp — > VVjj + X are presented. 
A program in Mathematica based on the structure of FeynCalc which automatically 
simplifies a set of amplitudes up to the hexagon level of rank 5 has been created 
for this purpose. We focus on two different topologies. The first involves all the 
virtual contributions needed for quadruple electroweak vector boson production, i.e. 
pp — > WW + X. In the second, the remaining "bosonic" corrections to electroweak 
triple vector boson production with an additional jet (pp — > VWj + X) are com- 
puted. We show the factorization formula of the infrared divergences of the bosonic 
contributions for WW and VWj production with V G (W, Z, 7). Stability issues as- 
sociated with the evaluation of the hexagons up to rank 5 are studied. The CPU time 
of the FORTRAN subroutines rounds the 2 milliseconds and seems to be competitive 
with other more sophisticated methods. Additionally, in Appendix A the master equa- 
tions to obtain the tensor coefficients up to the hexagon level in the external momenta 
convention are presented including the ones needed for small Gram determinants. 
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1 Introduction 



The incoming data from the CERN Large Hadron Collider (LHC) will require precise the- 
oretical predictions for a variety of signal and background processes. Processes with two 
vector bosons are of vital importance since they constitute background signals to Higgs and 
top physics as well as to physics Beyond the Standard Model (BSM) (for an overview see 
e.g. Ref. [I]). A tremendous effort has been carried out by the scientific community to com- 
pute these processes accurately. Next-to-Leading Order (NLO) QCD diboson production 
was computed in Refs. [2HI], and with an additional jet in Refs. [oTTTT]. Electroweak (EW) 
diboson production in association with two jets was computed in Refs. [T2HT1] . Recently, 
the NLO QCD corrections to W + W + jj + X [15J and W + W~jj + X [16] production were 
computed using the generalized D-dimensional unitarity framework for the calculation of 
the one loop virtual amplitudes. An impressive progress in the calculation of multi-parton 
one-loop amplitudes has been achieved not only applying new techniques [TT1 - HT] . but also 
based on traditional methods [IHrd]. In this paper, relying on traditional techniques, we 
compute some of the virtual corrections needed for pp — > VVjj + X production. We focus 
on some of the topologies that appear, namely QCD "bosonic" one-loop corrections to the 
diagrams contributing to pp —> WW + X and pp —> VWj + X production. The strategy 
followed is to collect Feynman diagrams with a given topology in groups which can be easily 
checked and reused in other processes. Thus, our aim is not only to provide results for 
VVjj + X production, including a second calculation of W + W~jj + X production using 
a different method, but also to provide a set of routines that can be used for many other 
interesting processes, such as W^j + X at NLO QCD [75]. It is known that in the calcula- 
tion of one-loop multi-leg amplitudes, the presence of small Gram and Cayley determinants 
might yield unstable results. In this paper, we show that the use of higher precision in the 
numerical determination of the tensor integrals (similarly to Ref. [76]) together with the use 
of Ward Identities to identify the unstable points, and to a minor extent, following Ref. [57] . 
the use of special routines for the determination of small Gram determinants for the C and 
D functions, solve this problem. 

The paper is organized as follows; in Section 121 the method used to perform the calculation, 
the different contributions computed and the tests performed to guarantee the correctness of 
the results are described. In Section [3j a discussion of the instabilities and the timing of the 
amplitudes computed are shown. The conclusions are given in Section |H In Appendix A, 
the tensor reduction routines used are presented. In Appendix B, numbers for the contribu- 
tions involving the hexagons are given, including proofs of the factorization of the infrared 
divergences and other additional numerical tests. Finally, in Appendix C, we give the color 
factors needed for testing the amplitudes. 

2 Calculational details 

In the calculation of the NLO QCD virtual corrections to pp — > VVjj + X, different sub- 
processes contribute; processes involving two quark pairs, e.g., uu — > ddW + W + and processes 
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with one quark pair, e.g., mu — > W + W gg 0. The latter can be separated into "one loop 
fermion" corrections and "bosonic" corrections, Fig. [TJ Among the "bosonic" contributions, 




Figure 1: "One loop fermion" contributions (left) and "bosonic" contributions (right). 

different topologies appear, Fig. |2j We will focus on diagrams with the first two topologies, 
i.e., one loop QED-like corrections and diagrams formed with one triple gluon vertex. The 
remaining ones will be discussed in future publications. 




Vi V 2 V 3 V 4 Vi V 2 V 3 Vi V 2 Vi V 2 



Figure 2: Different "bosonic" topologies. Vi stands for vector bosons emitted from the quark 
line. 

To compute these amplitudes, a function in Mathematica [77J using FeynCalc [78] which 
automatically simplifies a set of amplitudes up to Hexagons of rank 5 has been created. 
Throughout the calculation, the quarks are considered to be massless, the anticommuting 
prescription for 75 is used and we work in the Feynman Gauge. The result is given in 
terms of tensor integrals following the Passarino-Veltman convention [38], Appendix |Aj and 
written automatically in to FORTRAN routines. Nevertheless, the amplitudes can be eval- 
uated also in Mathematica with unlimited precision which is used for testing purposes. To 
achieve that, the scalar integrals, the tensor reduction formalism to extract the tensor co- 
efficient integrals, and also the helicity method described in Refs. [THIED] to compute the 
spinor products describing the quark lines of Figure |2] have been implemented at the FOR- 
TRAN and Mathematica level. For the determination of the tensor integrals up to the 
box level, we have implemented the Passarino-Veltman tensor reduction formalism [48J. We 
have also applied the LU decomposition method to avoid the explicit calculation of inverse 

1 Note that the NLO QCD corrections for pp — > W + W + jj + X production involve only sub-processes 
with two quark pairs. 
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Gram matrices by solving a system of linear equations which is a more stable procedure 
close to singular points. Finally, for singular Gram determinants, special tensor reduction 
routines following Ref. [63] have been implemented, however, the external momenta con- 
vention (Passarino-like) was used. The impact of these methods is discussed in detail in 
Section [31 For pentagons, in addition to the Passarino-Veltman formalism, the method pro- 
posed by Denner and Dittmaier [63], applied also to hexagons, has been implemented. For 
that, the recursion relations of Ref. [63] in terms of the Passarino-Veltman external momenta 
convention have been re-derived. This last method is used for the numerical implementa- 
tion at the FORTRAN level in Section [3J All the recursion relations used can be found in 
Appendix IA1 

In the following, a more detailed description of the method used is given. To extract the 
rational terms, two simplifications are done at the Mathematica level. To illustrate them, 
the amplitude of a simple vertex diagram is used, 




J (27r) d £ 2 {I+P1+P2Y {l+Pir 



<tP2 
fJ>2 <? 

First, a simple Dirac re-ordering manipulation is applied to explicitly contract repeated 
indices and to obtain all the terms proportional to the space dimension, D, coming from 
7 M 7 M = D contractions, i.e., 

D 

r f dH -( ^ aff^ _l„ ^ wIa.* \ 1 \ f d d l ^(p 3 )(/)7 M2 (/)7^^(Pi) , 



(2ir) d (l+ Pl +p 2 ) 2 (l+ Pl ) 2 l 2 



II 

(2.1) 

Second, all Dirac's structure containing the loop momenta is pulled to the right, such that 
resulting terms like // = I 2 are canceled against one of the integral denominators, 

II = Dv(p 3 )rMPi) [ ,tL mi ^ JtL I 

J {2TT) d l 2 (l+p 1 ) 2 (l+p 1 +p 2 ) 2 

oc Dv(p 3 )^u( Pl )I 2 (p 2 ) + ..., (2.2) 

where I 2 (p 2 ) = B (p 2 ) is the scalar two point function defined correspondingly to Eq. (1A.1I) . 
In this way we avoid possible additional terms proportional to D coming from the Lorentz 
structure of the tensor integral. More concretely, the same term II could result in terms 
proportional to D 2 : 



II = -Dv(p^ Ul ^ 2 ^ V2 u(p 1 



d d l l v H 



U 2 



(2ir) d P(i + Pl y(i + Pl + p 2 y 

oc —D v(p 3 ) luil ^ ll/2 u(Pi) x ifflfa&dg™ + ■■■ 

= -(2-D)Dv(p 2 )^u{p 1 )lS(p 1 ,p 2 ) + ..., (2.3) 
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with Iq (pi,P2) = Coo(pi>P2), a three point tensor coefficient integral of the tensor integral 
I3 111 * '(pi, P2) > defined by Eq. (jA.6j) . and obtained recursively from Eq. flA.8j) . Further sim- 
plifications are obtained using the Dirac equation of motion and rewriting the pair (/ • p^ 
as a difference of two propagators, e.g., (7 • pi) = (I + pi) 2 — I 2 , which are canceled against 
the denominators. The remaining rational terms stem from ultraviolet divergent tensor 
coefficients, which are treated independently within the tensor reduction routines. As an ex- 
ample, the finite contribution of the ultraviolet divergent tensor coefficient Coo, for massless 
propagators, following Eq. (1A.8|) is obtained by, 



a 



n=l 



,(UV) 



fin * V n=i / 

(2.4) 

where we have taken D = 4 — 2 e and series expand in e all the scalar and tensor coefficient 
integrals, e.g., 

B = B^ + \ B r\ (2.5) 

with B^ v ^ = 1. After these steps, the amplitudes, M. v , in terms of scalar and tensor 
coefficient integrals can be written by, 

M v = M° =A + (D-4)M° R , (2.6) 

where Ai^ =A is the amplitude that one would obtain performing the Dirac algebra ma- 
nipulation in four dimensions, D = 4, and M.® R contains the rational terms and vanish 
in Dimensional Reduction (DR). To get this expansion, we have only rewritten the space 
dimension D as D = d + 4 such that for example Eq. (12. 2 ft is given by, 

Dv{p 3 )Y 2 u{Pi)B {P2) = 4v{p 3 )YMPi)Bo{P2) +d v{p 3 )YMPi)Bo(P2) • (2.7) 

V V ' S v ' 

Both M.® =A and M.® R are decomposed in terms of, 

M {d=a,dr) = SM iiT Flj, (2.8) 

where SM^ T is a basis of Standard Matrix elements corresponding to spinor products de- 
scribing the quark line of Figure |2] which are computed following the helicity method [791180] 
with a defined helicity, r. Flj are complex functions which are further decomposed into 
dependent and independent loop integral parts, 

Flj = ^F fc T fe (e(p n ) -p^eipi) -e( Pi )). (2.9) 
k 

Tfc is a monomial function at most for each polarization vector e{pS). contains kinematic 
variables (p, -Pj), the scalar integrals (B = I2,C = h, D = J 4 )q and the tensor integral 

2 Note that throughout the paper to name the different n-point function integrals, /„, the alphabetically- 
ordered labeling frequently found in the literature is also used. 
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coefficients (Bij,Cij,Dij,Eij,Fij). The latter obtained numerically from the scalar integral 
basis following the recursion relations of Appendix A. To illustrate this notation, the following 
example is given, 



Fli 



~>(P4)p 2 P+u(pi) C (px,p 2 ) e(pa) • e(pa) = v(p 4 )p 2 P+u(pi) C (p 1 ,p 2 ) e(p 3 ) • e(p 2 ) = SM^+Fh 



SMi + T i 

(2.10) 



P + is the positive helicity projector, P + = We note here that this parametrization 

is convenient for example to evaluate the amplitude for different polarization vectors or for 
performing gauge tests since some of the functions will remain unchanged. Finally, the full 
result is obtained from A4® =i and At^ R using the finite and the coefficients of the l/e n 
poles of the scalar and tensor coefficient integrals (see e.g. Eq. (12.51) ). From Ai® =A , we get 



Ml At 

e 2 



M^ = M v + ^ + (2.11) 



where Ai v is the finite contribution obtained using the finite pieces of the scalar and tensor 
coefficient integrals including the finite contributions from rational terms arising in ultraviolet 
tensor coefficient integrals, Eqs. (12 .4p and flA.8[) . Ail and Ail are obtained from the 1/e and 
1/e 2 pole contributions, respectively. Similarly, from (D — A)Ai DR , one gets, 

(D-4)Ai° R = Af v + -^, (2.12) 

where Af v and Af* are obtained from the 1/e and 1/e 2 poles, respectively. Indeed, Af* is 
analytically zero since rational terms in one-loop QCD amplitudes, omitting wave renormal- 
ization graphs (WRF) are of UV origin [ST]. Eq. ( 12 . 1 2[) will allow us to check this statement 
numerically. Additionally, up to the pentagon level, the divergent part is computed in Di- 
mensional Regularization analytically in Mathematica. A library containing all the divergent 
tensor coefficients for massless particles including pentagons of rank 4 has been created with 
this purpose. This allows us to obtain also the terms of Eq. (12. 12j) and the divergent parts 
of Eq. (12. lip analytically up to this level. 

It is known that the IR divergences depend on the kinematics of the external particles 
involved in the process, i.e., whether they are massless/massive on-shell/off-shell particles. 
Nevertheless, since we reconstruct the divergences numerically, to compute generally the 
amplitude of a particular diagram for a given helicity, r, it is sufficient to consider off- 
shell vector bosons, Vi in Figure |2j This basically means that all pf terms are kept and the 
transversality property of the vector on-shell bosons (e(p x )-p x = 0) is not applied analytically, 
allowing us, the latter, to perform gauge tests (e(p x ) — > p x ) for on-shell massive particles. 
We note here that the polarization vectors can be understood as generic effective currents, 
t^iPx) J^ffiPx), which can include the leptonic decay of the vector bosons or physics BSM. 
Thus, denoted by Aiv 1 v 2 v 3 V4 the amplitude of the first diagram of Figure |2] would be given 
by, 

M*w«,r = ^W^^Cf^^f, (2.13) 
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where go is the strong unrenormalized coupling, Ai l J is the corresponding diagram with 
off-shell vector bosons (or effective currents) with color indices ij and given in terms of 
Eq. f)2.6p which is computed only once and for all vector bosons, V { e (7, W, Z,g). C^ 2 ^ 4 
is a color diagram dependent factor, e.g., C^ gi = (T a )ij(Cp — 1/2Ca) and from now on the 
color sub-indices will be omitted, g™* is the coupling following the notation of Ref. [80] with 
e.g., gl f = \e\Q f , g_ f = \e\/(V2sm9 w ), gj = |e|(T 3/ - Q f sin 2 9 w )/(sm9 w cos9 w ), where 
9w is the weak mixing angle, T 3 j is the third component of the isospin of the (left-handed) 
fermions, and |e| the electric charge. 

This procedure is advantageous for several reasons. First, when considering off-shell vector 
bosons, for example, diagrams with the same QED-like topology as the first one of Fig. [2] 
will give us, already, all the corrections needed for EW quadruple vector boson production, 
pp — > WW + X and part of the pp — > VWj + X and pp — » Wjj + X contributions, 
independently whether we are considering massless or massive particles or the bosons emitted 
from the quark line are gluons. Cjj V2V3V4 will associate the correct color factor to each 
diagram depending whether one is considering pp — > WW + X, pp — » VWj + X or 
pp Wjj + X and g^f the corresponding couplings. Also, the leptonic decay of the vector 
bosons and new physic effects can be immediately incorporated by using the appropriate 
effective currents. Second, cross term related diagrams are obtained from the same analytical 
amplitude just by permuting the vector boson momenta and polarization, and thus, e.g., 
from the 24 cross related diagrams to the first one of Fig. [2j only the permutation depicted 
has to be computed. Moreover, additional checks can be implemented, as for example the 
factorization of the IR divergences against the born amplitude for different processes starting 
from the same analytical structure. In the following, we will refer to the first topology of 
Fig. [2] as contributions to pp — > WW + X and, to the second as "bosonic" contributions 
to pp — > VWj + X since they appear for the first time in these processes. 



2.1 Contributions to pp WW + X 

In this section, all the "bosonic" QED-like loop corrections along a quark line needed for 
pp — > Wjj + X production are considered. These contributions can be classified by I) 
virtual corrections along a quark line with two vector bosons attached (first diagram of 
Figure [3]) II) virtual corrections with three vector bosons attached to the quark line (second 
and third diagrams of Figure [3]) and III) virtual corrections along a quark line with four 
vector bosons attached (last diagram of Figure [3]). The strategy followed is to search for 
a minimal set of universal building blocks from which every amplitude can be constructed. 
We use the effective current approach described and applied in Refs. [7] [T0l[TT| I52TI81)] . As 
illustration, the first diagram of Fig. [3] is used. This can be written as, 

^ViV 2 v 3 V4,t = JvfJvz-Mp^T = Mv*v*,r, (2-14) 

where the color indices have been omitted. Here, Jyl and Jyl represent effective polarization 
vectors including finite width effects in the scheme of Refs. [T2"ll8~7] and propagator factors, 
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V 3 V 2 






Vn V\ V?> V 4 



v l v 2 v 3 v 4 



Figure 3: "bosonic" one loop QED-like topologies appearing in the calculation of the virtual 
contributions for pp — >• VVjj + X production, with Ve (W ± ,Z,7). 



e.g. 



Ml 



Ml* 



IM V *T V * 



IM V *T V * 



1 V^ViVs 



(2.15) 



with r^*, the width of the V* vector boson, and Ty, ViV3 , the triple vertex which can also 
contain the leptonic decay of the EW vector bosons including all off-shell effects or BSM 
physics. In this manner, we can treat diagrams with triple vertexes as external off-shell 
legs (in the sense that they will have the same IR behavior). We can then concentrate 
in computing, instead of A Vi v 2 v 3 V4,t, the virtual correction to two massive vector bosons 
attached to the quark line, M.y^y*^-, or equivalently A^ MlM2)T , where the polarization vectors 
or effective currents have been factored out. Thus, the strategy will be to combine in groups 
all the virtual corrections to a given born-amplitude configuration, independently of the 
effective currents or polarization vectors attached to the quark line. Furthermore, the order 
of the gauge bosons are fixed and the full amplitude will be recovered by summing over 
the physically- allowed permutations. Thus, three universal virtual contributions are left, 
corrections to a born amplitude with two, three and four vector bosons attached to the 
quark line for a given order permutation. 

I) The virtual corrections to the Feynman graph with two vector bosons V\ and V 2 attached to 
the quark line for a given order-permutation with incoming momenta k\ and k 2 are depicted 
in Figure H] with kinematics, 



?(pi) + 5(p 2 ) + V^h) + V 2 (k 2 ) -> 0. 
The sum of the four amplitudes with a given helicity is written in terms of, 

i 

(4tt) 



(2.16) 



n 2 

M Vl V 2 ,r ~ 9r J 9r TJZ^Z^ C (n) M 



(n) 

ViV 2 ,r 1 



71=1 



M 



(n) 

VlV2,T 



A/f» 4- Kf 



r(n) 



J^WV^t y> Jvl ViV2,r 



(2.17) 



where A4viV 2 ,t is called from now on "boxline" contribution. Although, we are interested on 
off-shell vector bosons emitted from the quark line (effective polarization vectors), the diver- 
gent contributions and A/y"} 2 T were also computed analytically for the different kinematic 
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Vi v 2 V, v 2 v l v 2 

(2) (3) (4) 

Figure 4: Virtual corrections for a fermion line with two vector bosons attached, V\(ki) and 
V 2 {k 2 ) in a given permutation. The sum of these graphs defines M.ViVi,r m Eq. (12.171) . 

configurations. These configurations are obtained by considering the vector bosons massless 
or massive and all possible ordering permutations, i.e., {k\ — 0, k\ — 0), (kf = M\,k\ = 
0), (kf = 0, k\ = M 2 ), {k\ = Mi, k\ = M 2 ). From the analytical calculation, it is confirmed 
that the rational terms arise from diagrams with ultraviolet divergences, thus, Nyy T is 
zero. For configurations with massless vector bosons, kf = 0, the transversality property of 
the bosons must be used, ki ■ e{ki) = 0. A general proof of this statement was presented in 
Ref. [SI] . Moreover, Myy^ T is zero since there are not 1/e 2 UV poles. For the self-energy 
diagram, M.yy 2 T =0 since the divergences are only 1/e. Finally, for M.yy 2 T , there are not 
collinear and soft singularities simultaneously given rise to 1/e 2 poles, so that, M v % ^=0. 

II) The virtual corrections to the Feynman graph with three vector bosons V\, V 2 and V 3 
attached to the quark line for a given permutation are depicted in Figure |5] with kinematics, 

?(pi) + <i(P2) + Vi(h) + v 2 (k 2 ) + v 3 {k 3 ) 0. (2.18) 

The sum of the eight diagrams with a given helicity is written in terms of, 

2 8 

KA — r,Vlf „V2f n V 3f S ° V^lWs M<") 

M Vl VzV Z ,T ~ 9r 9r 9r 7T32 2_^ L » M ViWs,T I 

^ ' n=l 

M-vmv^t ~ ■ M v 1 v 2 v 3 ,T+ J ^v 1 v 2 v^r + + 2^ a ' ^ IA ^> 

where M.ViV2V 3 ,t is called from now on "penline" contribution. We have computed also the 
divergent contributions and Nyy 2 y 3 T analytically for the eight different kinematic configu- 
rations, (kj = 0, k\ = 0, kl = 0)) \k\ = Mi, k\ = 0, k\ = 0), (kj = M u k\ = M 2 , k\ = 0), . . .. 
From the analytical calculation, it is verified that once the transversality property for mass- 
less on-shell particles is applied, the rational terms arise from diagrams with ultraviolet 
divergences. Therefore, Nyy^y 3 T and My^y T are zero. Analogously to the boxline, for 

self-energies, My V fy 3 T =0 as well as My V ~yl T =0. 

III) The virtual corrections to the Feynman graph with four vector bosons V\, V 2 , V 3 and V4 
attached to the quark line for a given permutation are depicted in Figure [6] with kinematics, 

q(pi) + <i(p2) + v^h) + v 2 (k 2 ) + v 3 (k 3 ) + v 4 (h) o. (2.20) 
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V 1 V 2 V 3 Vi v 2 v 3 

(7) (8) 

Figure 5: Virtual corrections for a fermion line with three vector bosons attached in a given 
permutation. The sum of these graphs defines M.ViV 2 v 3 ,t in Eq. (12.191) . 

The sum of the thirteen diagrams with a given helicity is written in terms of, 

2 13 

M vvvv - a Vif a Vif a Vaf a V*f 9o p V 1 V 2 V 3 V4 A J W 

^iViV 2 v 3 Vi,T — y T y T Vt y T r± K y / jin\ j 1 ViV 2 v 3 v 4 ,t > 

^ ' n=l 

A/fW - jCfW 4- I lljjjW^I I • /V VlV r 2V r 3V r 4,T / 991 N 

1=1 

M-V\V 2 v 3 Vi,T w ih be called "hexline" contribution in the following. For the hexline contri- 
bution, the divergences are not computed analytically. The explicit and naive recursive 
reduction and simplification of the tensor integrals for hexagons in terms of 1/e poles ex- 
ceeds the capacity of Mathematica with 4 GB of memory RAM. Therefore, Eqs. ( I2.11|2.12p 
are used to obtain numerically the different contributions both at the Mathematica and 
FORTRAN level. We have checked with up to 30000 digits of precision in Mathematica 
for the 16 different kinematic configurations, i.e., (k\ — 0, k\ = 0, k% = 0, k\ — 0), {k\ = 
Mi, k\ — 0, k\ — 0, k\ = 0), (k\ = Mi, k\ = M 2 , k\ — 0, k\ — 0), . . . and several phase space 
points that ^Vivtv 3 Vi r * s zero ^ or an possible kinematic configurations (for massless particles 
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Vx V 2 V 3 V 4 Vi V 2 V 3 V 4 Vi v 2 v 3 v 4 

(11) (12) (13) 

Figure 6: Virtual corrections for a fermion line with four vector bosons attached for a given 
order permutation. The sum of these graphs defines M.y^VzV^r in Eq. (12.2 ip . 

the transversality property must be used) as well as •N'viW/zVa t = ■^Vi^yfvkr = ^- At the 
FORTRAN level for non-singular points, the proof works at the working precision level. 



2.1.1 Checks 



In this section, the tests performed in order to ensure the correctness of the calculation of 
the 1 loop diagrams/contributions are explained. We consider electroweak vector boson pro- 
duction (Vi,V2, ^3; Vi) ^ (W^Z^). Then, the color factor for all diagrams and contributions 
is proportional to Cp. This is enough to test not only the individual diagrams but also the 



11 



different contributions. For electroweak vector bosons, the divergent terms for the boxline 
and penline are known analytically, the general form generalizes for n vector boson emission 
and reads, 



M 



Vl...Vn,T 



n Vif n v n f r 



47T 

r(i + e) 



M v + 



10 - D 

2c 



(2.22) 

(^)e(W ± ,z, 7 ) 



where .My VnT is the corresponding Born amplitude with helicity r and s the square of 
the partonic center-of-mass energy. If we use cr(— s) in Eq. flA.2[) for the definition of the 
scalar and tensor integrals and D = 4 — 2e, then, the pre-factor (^) e T(l + e) of Eq. (12. 22ft 
is reproduced and each of the terms of Eq. f !2.22|) can be identified, for example, for the 
hexline contribution, Eq. (12. 21 j) . with the different finite and pole terms as, 



13 



13 



13 



n=l 



(n) 

VlV 2 V 3 V4,T 



M v , J2 M 



n=l 



l,(n) 

ViV 2 V 3 V4,t 
13 



— 3 KA B 



J2 M 

n=l 



2,(n) 

VlV 2 V 3 V4,T 



2-A / ly 1 y 2 y 3 y 4 , T , 



(n) 
V a V 3 Vi,T 



-M B VA/ 1,(n) 



(2 



n=l 



n=l 



These relations are very important for testing purposes. Note that the left hand quanti- 
ties for each of the above lines are obtained numerically from the same analytical expres- 
sion, Eqs. (I2.11f2.12p . They only differ in which terms of the l/e n expansion of the scalar 
and tensor integrals are used. Therefore, the numerical check of the factorization of the 
singularities and rational terms provides a strong check of the correctness of the finite terms 
■MviV 2 v 3 v 4 t m We have checked the factorization formula analytically for the boxline and pen- 
line contributions for all different kinematic configurations. To cast the singularities in this 
form for massless on-shell vector bosons, (7), the transversality property of the particle must 
be used. Otherwise, additional singularities proportional to the product k { ■ e{kj) appear, 
with ki the on-shell massless particle momentum. Thus, we have checked that the presence 
of additional massless particles does not introduce new singularities or new rational terms 
in agreement with Ref. [81]. For the hexline contribution, we have checked numerically with 
high precision in Mathematica the factorization formula for the 16 different kinematic con- 
figurations and for different phase space points . The coefficients multiplying the poles of 
Eq. (12.221) are obtained with up to 30000 digits of precision at least (at the FORTRAN level, 
for non-singular points, the proof works at the working precision level). Nevertheless, we 
will assume the result to be analytic due to the high precision achieved. An analytical proof 
can be obtained using the method described in Ref. [SB] which is not automatized within 
our approach. 

Concerning the origin of the rational terms in Eq. (12.221) . we note that for electroweak 
boson production, the wave renormalization functions (WRF), which are zero in Dimensional 
Regularization for massless on-shell quarks, together with A4v 1 ...v„,t are UV finite. Therefore, 
all the divergences in Eq. (12. 22ft . after adding the WRF, become of IR origin including the 
terms containing the rational factor, D/(2e), of UV origin. This is clear since the WRF 
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are zero due to the cancellation of IR and UV poles. Treating the UV and IR divergences 
separately, one observes that the UV poles of the WRF cancel the rational terms of Eq. (I2.22p 
and only IR divergences remain (see also Ref. [ST]). 

The factorization proof already provides an important check of the correctness of the calcu- 
lation. We can use the factorization formula to perform an additional test. The factorization 
of (— s)~ e from the scalar integrals introduces an additional dependence in Ai v on this vari- 
able. If we consider (-s) as an independent scale energy variable, — s = fiQ, and series expand 
(yU )~ e m e in the second piece of Eq. ( I2.22p . we obtain a new finite term, Ai' v , of the form[f|, 

M' v = Mv(^o) + /(#>) • Mb, with /(// ) = - (log 2 (/i ) - 31og(/x )) , (2.24) 

such that the log(/xo) terms compensate the no dependence propagated in all the scalar and 
tensor coefficient integrals in the complex _A/fy(/U ). Thus, A4' v is /i independent. This 
fact was satisfied for non-singular points in the three contributions at the working precision 
level in the FORTRAN code (for the hexline see Tab. [T2|) and at least with 30000 digits in 
the Mathematica code. Moreover, we have implemented Ward identity tests for the virtual 
corrections in FORTRAN and Mathematica at different levels of complexity: 

1) At the level of single diagrams. 

2) For the MviV 2 ,t, Mv^v^t, M-ViV 3 v 3 V4,t contributions. 

3) At the level of gauge invariant quantities. 

4) Subset of amplitudes invariant for a specific replacement (e(pfc) = — > Pk)- 

5) Specific contractions that make the contributions to vanish. 




Figure 7: J> 2M3M4M5 (pi,p 2 ,P3,P4,P5) and S^ 2tl3tl4 (p 1 ,p 2 ,p 3 ,pi) of Eq. (|2.25|) and Eq. (j2.26[). 

3 Note that A4v(no) and My, which is A4v(no — lGeV), have the same analytical structure. They only 
differ in the input scalar integrals. 
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1) Under the replacement of the polarization vector by its momentum (e(pk) = tk — > Pk), 
relations among the different diagrams/topologies can be found. As illustration, one can 
consider the hexagon of Fig. 

v < \ f <]dq 1 " 1 1 1 1 1 

(2.25) 

where py = Ylk=iPk- Contracting one of the open indices by the corresponding momentum 
and expressing the contracted gamma matrix as the difference of two adjacent fermionic 
propagators, the hexagon is reduced to a difference of two pentagon integrals, 

P2,P3,P4,Ps) = £i*miMi(pi,P2+Pa,PA,P5) -£w**«i(Pl+P2>P3,P4,P5), 

P2,P3,P4,P5) (Pi, P2, P3 + P4, Ps) - £ M 2M4M5 (Pl> P2 + P3, P4, Ps), 

Pi A F^z^{Pl,P2,P3,P4,P<o) = £/J,2fi3Hs(Pl,P2,P3,PA+Ps) ~ £ M 2M3Ms (Pi, P2, P3 + P4, Ps), 

Pf U2V3V-Altb(P^lP2lP3lP4:-,Pb) — £/12M3/m(P1> P2,P3,Pa) - S llalMllH (pi,p2,P3,P4 + Pb)- (2.26) 

where S^^^ represents the pentagon diagram of Figured defined similarly as Eq. (12.251) . 

2) For the contributions, assuming EW vector boson production, equivalent relations can 
be obtained. First, we factor out the couplings and the polarization vectors, 

M Vl ...V n ,r(Pu ...,Pn)= 9r lf ■ ■ ■ 9r nf tfM) ■ ■ ■ <^(P»)^U...A*,r(Pl, • • ■ , Pn) ■ (2.27) 

where ■M At2 ... A t n , r (Pi> • • • >Pn) represents the un-contracted contributions and we have explicitly 
written down the momentum's dependence following the convention of Figure [7J Note that 
the color factors, C^'" Vn , are still present in A^ M2 ... MnjT (pi, • • • ,p n )- Then, Eq. (12.261) holds 
for the hexline and penline contributions with the replacements, 

^mnjni M-mnjni^- (2.28) 

Analogously, the penline contribution can be checked against the boxline contribution ob- 
taining similar relations. Additionally, to check the boxline contribution, we computed the 
quark line with one vector boson attached, which we called vertline, Ai^ 3lT , and is formed 
by a simple vertex graph. These relations are satisfied independently both for the finite and 
the divergent parts of each contribution. 



3) Vertline is gauge invariant under the replacement of e Mi (pj) — > pf* which means that the 
equivalent form of Eq. (12.261) gives zero. To build gauge invariant quantities for the other 
un-contracted contributions and for EW vector boson production, the cross diagrams must 
be considered (See Fig. E as illustration). Thus, the sum of direct and cross terms under 
the replacement e Mi (pj) — > pf' is gauge invariant: 

Pi' (M ll2fl3iT (p l ,p 2 ,P3) + M ll3ll ^ T (p 1 ,p 3 ,p 2 )) = 0, 

Pf (- A/ ^M2M3M4,r(Pl,P2,P3,P4) + M^^Pl, P3, P2, Pi) + 4 Other) = 0, 

Pf (- A/ ^M2M3M4M5,r(Pl,P2,P3,P4,P5) + A / f M 3 M2At4M5iT (pi,P3,P2,P4,P5) +22 Other) = 0. 

(2.29) 

4 Note that for processes with W vector bosons not all the permutations are physically allowed. 
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Figure 8: Direct and cross term diagrams which correspond to permuting the momenta and 
the polarization vectors in the boxline contributions yielding a gauge invariant subset for 
Vie {i.,Z). 

4) For the un-contracted penline and hexline and considering EW vector boson production, 
there are smaller subsets which give zero ("gauge invariant") for a given replacement {e{pk) = 
£fc — >Pk)- For the un-contracted penline A^ At2M3/i4jr (pi,P2,P3,P4), a "gauge invariant" subset 
for a specific replacement, e.g., e(j>2Y 2 — > P2, 2 is obtained by permuting the position of the P2 
momentum and the corresponding Lorentz index and keeping the relative order of the other 
vector bosons fixed, i.e., 

{■M- fl 2 tJ,3fl4,T{Pl, P2, P3, Pi) + -Mm3M2M4,t (Pl,P3,P4,P2)) = 0. 

(2.30) 

The same can be proved performing the corresponding replacement (e(pk) = e& — > Pk) for 
the contracted contributions which include the couplings g% n * , subject that V n G (7, Z). 

5) In addition, we have checked analytically and later on numerically that there are some 
combinations of replacements for the polarization vector which make for EW vector boson 
production the contributions to vanish, e.g., 

(Pf (P2+P3) M3 ||(P2+P3)^f) (AW 3! t(P1,P2,P3)) = 0, 

P2 2 (P2+P3+P4Y S pT (M fl2flsfl4! r(pi,P2,P3,Pi)) = 0, 

(P2 2 (P2 + PsT 3 (P2 +P3+P4+ PSY 3 PT I W (P2 +P3+P4+ PsY* (?4 + PsY'Ps) 

(M II2fl3II4fl5 , T (p 1 ,P2,P3,P4,P5)) = 0. 

(2.31) 

These tests have been implemented at the FORTRAN and Mathematica level. At the 
Mathematica level a precision of 30000 digits was achieved. At the FORTRAN level, for 
non-singular phase space points, these tests are satisfied at the working precision level. The 
finite boxline and penline contribution have essentially the same analytic expressions found 
in the calculation of the NLO QCD corrections in vector boson fusion processes, qq — > Vqq 
and qq — > VVqq, discussed in Refs. [12] and [15] . respectively. We have checked that our 
results agree at the double precision level. The use of modular structure routines, as the 
above presented, has been proved to be an advantage in the program VBFNLO [89] since 
once a structure is computed and checked it can be reused for different processes. For 
example, the boxline and penline contributions computed here have been used to compute 
the NLQ QCD corrections to W ± W ± Z [82], W + W~-f (ZZ-y) [B3J, W ± Z 1 [M], [85], 
W ± jj 0EEI], W ± Zj pUiEE] and also Hjjj [90] production. They are publicly available as 
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part of the VBFNLO package together with the tensor reduction routines up to the pentagon 
level, excluding the routines for small Gram determinants. 




V 1 V 2 



Vx V 2 



Vi v 2 v 3 



Vi v 3 v 2 



Figure 9: Topologies appearing in the calculation of the virtual contributions for pp — > 
VVjj + X production, with V<G (W ± ,Z,7) . 



2.2 "Bosonic" contributions to pp — > VVVj + X 

One loop non-abelian corrections involving one triple gluon vertex are considered in this sec- 
tion, Fig. El We follow the same strategy used in the previous section and classify the virtual 
corrections to three different groups depending whether they contribute to I) qqV*g* — > 0, 
II) qqViV 2 g* — > or qqV*V 3 V4 — > or III) qqViV 2 V^V A — >■ 0. In analogy to the previous 
section, for the last two, only diagrams with a specific permutation of the vector bosons are 
considered. Generally, the momentum square for the off-shell legs V{ and g*, indicated with 
a star, is not zero. This is important since it results in different IR divergences for individual 
graphs. From now on, to avoid confusion, the star is only kept for the gluon, but we remind 
the reader that at the programing level all vector bosons are considered to be off-shell as 
described in Section [2] since the divergences are reproduced numerically, Eqs. (I2.lip2.12p . 

I) Non-abelian corrections to the Feynman graphs with one vector boson V\ attached to the 
quark line are depicted in Figure [10] with kinematics, 



q(pi) + q(p2) + V x {k x )+g*{k 2 )^U. 



(2.32) 



The sum of the three graphs with a given helicity is written in terms of, 




2 



M 



M 




(2.33) 



i=i 



where M.y T is called from now on "boxlineNoAbe" contribution. We have computed the 
divergent contributions and J\f Vi analytically for the 4 different kinematic configurations. 
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o 
o 
o 
o 




Figure 10: Virtual corrections for a fermion line with one vector boson attached, Vi(ki) 
The sum of these graphs defines M.y T in Eq. (12.331) . 



From the analytical calculation using the transversality property for massless particles, it is 
confirmed that My and Mi 



'Vl,T 



are zero. 



II) The virtual corrections to the Feynman graph with two vector bosons V\ and V 2 attached 
to the quark line for a given permutation are depicted in Figure [11] with kinematics, 



?(Pi) + fe) + Vrih) + V 2 {k 2 ) + g*{k 3 ) -> 0. 
The sum of the six graphs with a given helicity is written in terms of: 



(2.34) 



M 9 y lV2>T 



M 9 *' (n) 

JVl V 1 V 2 ,r 



V\f v 2 f 9o 
9r 9r 9o 



',ViV 2 KA9*,{n) 



U ™ " 'ill 



n=l 



J\sl9*,(n) KT9*,(n) , 



JV ViV 2 ,t 



M 



9*,i,(n) 

ViV 2 ,T 



(2.35) 



where M.y x y T is called "penlineNoAbe" contribution in the following. We have computed the 

divergent contributions and N Vl y™\ analytically for the 8 different kinematic configurations. 

It is verified that Nyyl"^ and Nyy^r are zero (^ ne transversality property for massless 
particles must be used). 



Ill) The virtual corrections to the Feynman graph with three vector bosons V\, V% and Vs 
attached to the quark line for a given permutation are depicted in Figure [12] with kinematics' 5 ' 



?(Pi) + ?(P2) + Vx{k x ) + V 2 {k 2 ) + V 3 (k 3 ) + g*(h) -)• 0. 



(2.36) 



5 Note that we do not need to consider an off-shell gluon in this case since for pp — > VVjj one of the Vi 
must be also a gluon in Eq. (j2.36[) . 
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9 

o 
o 
o 




v l v 2 

(3) 





9 




Vi 



Vi 



(5) 



Vl 



Figure 11: Virtual corrections for a fermion line with two vector bosons attached, Vi(ki) and 
^2(^2) for a given order permutation. The sum of these graphs defines A4.y lV2T in Eq. (12.351) . 

The sum of the ten diagrams with a given helicity is written in terms of: 



M 9 *' {n) 

jvv VJ/ 2 V 3 ,t 



2 10 

a Vif a V 2 f V 3 f \- C 9,ViV 2 V 3M g*,(n) 

^ > n=l 



M 



9*,{n) 
V 1 V 2 V 3 ,t 



V 1 V 2 V 3 ,t 



+ 



i=l 



2 \A9*,i,(n) 
Jvl VxV 2 V 3 ,n 



(2.37) 



V1V2V3 



r is called "hexlineNoAbe" contribution in the following. We have checked with 



high precision numerically in Mathematica that once the transversality property for massless 



particles is used, Af{ 



■g* (1-6) 

ViV 2 V 3 ,t 



and Afi 



ViV 2 V 3 ,t 



are zero for the 16 different kinematic configurations 
(at the FORTRAN level for non-singular points, this proof works at the working precision 
level) . 



18 




Eq. (E37D. 



2.2.1 Checks 

The factorization of the infrared divergences for one or two electroweak vector boson pro- 
duction in association with one jet, gu —> uV\ or gu — > uViVz, is known analytically [T|l9Tj . 
We will make use of this fact to perform powerful tests to our contributions. For one gluon 
emission, the color for all the non-abelian contributions is proportional to Ca, and two color 
factors proportional to Cp and to Cp—\j1C a (see Fig. [13]) appear in the abelian contributions 
computed in the previous section. The relevant color factors for this proof for all diagrams 
of each abelian contribution are given in Appendix O We will concentrate in the hexline 
contributions since the penline and the boxline contributions are checked similarly. To verify 
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oc Cf-| C A oc Cf~\ Ca oc Cp 

Figure 13: Color factors that appear when one gluon is emitted from the quark line. 

the factorization of the divergences against the born amplitude for (Vi, V 2 , V 3 ) G (W ± , Z,j) 
and for an on-shell gluon, it is enough to consider all the diagrams mixed under the QCD 
gauge group for a given order combination of the vector bosons, (Vi, V2, V3), 

■MgViV2V3,T + Mv ig V2V 3 ,T + Mv!V2gV 3 ,T + M-V\V 2 V 3 g,T + -M-VlVlVz,T- (2.38) 

The same is proved for the other six permutations (Vi V2 V3). This combination 
under the replacement e(/c 4 ) — > fc 4 , given the kinematics defined by Eq. (I2.36|) is "gauge 
invariant" . The virtual contributions factorize on the Born amplitude generalizing the result 
of diboson plus jet production. The sum of the amplitudes, including the counterterms (CT), 
in Dimensional Regularization are given by, 



MgV-iVzV^T + ■M.V 1 gV 2 V 3 ,T + M-V x V 2 gV z ,T + ■M.V 1 V 2 V 3 g,T + ■^V 1 V 2 V a ,T + CT 



47r v ~' \2 \\ —u J \ -t J J v e 2 e 
e 2 e 4e 



+ 2 ('^V ("% - * + ^^i>) I + (K„ K„ VJ) e (W*, Z, 7 ), (2.39) 



with 



CT = -M??^r{l + e){i*)& l q = 3 -C F , lg = ^C A - 2 -T R N f , (2.40) 

T fl = 1/2, C A = N, and C F = (N 2 - l)/(2iV) in SU(N) gauge theory. The number of flavors 
is Nf = 5, s is the square of the partonic center-of-mass energy, s = (pi + P2) 2 , w = (pi + /C4) 2 
and £ = (p2 + ^4) 2 , given the kinematics defined by Eq. (12.361) . A^f is the Born amplitude 
formed by four tree level graphs, the corresponding born amplitudes of each of the "hexline" 
contributions, M? = My x y 2 y zgT + M Vl y 2gV3T + M Vig y 2 y 3T + M^ Vl y 2 y 3T . Given our definition 

for the scalar and tensor integrals, Eq. (1A.2j) . M.y is defined through, 
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EfpVMVsg Tj(n) nV^gVs V7(n) nV ig V 2 V 3 T~7(n) r gV x V 2 V 3 7j{n) \ 

\^{n) JVi ViV 2 V-ig,T ~+~ u (n) JVl ViV 2 gV 3 ,T u (n) Jvl V ig V 2 V 3 ,r L » "^gViV&V^T J 

n=l 

10 

+ E C uJ lV2V3M v$ 2 v 3 ,T = M V + M B T f{u, s, t, /x 2 , C F , C A , 7 ,), (2.41) 
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where /(it, s, t, /i 2 , Cp, Ca, J g ) is given by the finite terms resulting of the epsilon expansion 
ofEq. (E^D, 

nu,s,ty,C F ,C A , 7 ,) = -C F (log - log + | (log g) + 

+ ^)) + ^K*&)-*G)-*(?)) + 

+ 2^(- 2 ^)+^) +1 °s (?)))■ < 2 - 42 ' 

From the analytical form of the factorization formula, we can find relations for the divergent 
and rational terms. For the terms proportional to Ca, we get 

— ( V M Un) + V M Un) + V M Un) + V M Un) 

2 I Z^ y Vl ViV 2 V 3 g,T + w ViV2 9 V 3 ,t + 1 Vi S V2V3,t + Z^ l sViV 2 V 3 ,r 

\n=l,3,6,10 n=l-3,5,6,10 n=l-5,8 71=1,2,4,7 y 

+ £ - (- >og fe) + log fe) + log fe) ) , 



Ji J V /•* / V A* 

— I V 7W 2 ' (n) + V 7W 2 ' (n) + V M Un) + V M Un) 

\n=l, 3,6,10 n=l-3,5,6,10 n=l-5,8 71=1,2,4,7 y 



10 

+ VA/f 9 ' 2,{i) - -A/f B 
+ Z^ J l viv 2 y 3 ,T _ yvi T ' 

-1 



+ 



-1 

IT 



+ 



Z^ J ^ViV 2 v 3 g,T + Z^ • /v Viy2gV 3 ,r + Z^ ■ / Vi<?V2V 3 ,t + Z^ ^aViv^r 

v n=l,3,6,10 n=l-3,5,6,10 n=l-5,8 71=1,2,4,7 / 

10 

E Art/, CO _ n 
• A ViV 2 V 3 ,t ~ U ' 

i=i 

V JV 1,(n) + V 7V 1,(n) + V Af Un) + V JV 1,(n) 
/ , J ^ViV 2 v 3 g,T + Z^ • /v Viy 2 gV3,T + Z^ • /v Vi 9 y 2 y 3 ,T + Z^ ^aViV^r 

,n=l,3,6,10 71=1-3,5,6,10 n=l-5,8 71=1,2,4,7 J 

10 

XXw = °. (2.43) 
i=i 
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and for the Cp terms 

13 



E 



(MlM^r + M\%U.r + Mtlv„, + M^ r ) = (-3 + 2 log (^)) 



71=1 

13 

/ j y vl ViV 2 V 3 g,T + Jvl ViV 2 gV 3 ,T + •^Vis^Vs.t + -^'gViV^r ^ zyw r > 

n=l 
13 



E 



(V (n) + A/" (n) 4- \f (n) + A/" (n) ^-_A/f B 



=1 
13 

E(a^<»> +JV 1 ' (n) + A/" 1,{n) + Af 1,(n) 1— (2 441 



We have checked numerically with high precision in Mathematica the factorization formula 
for the 8 different kinematic configurations of the hexlines once k\ = is fixed, i.e., (kf = 
0, fc| = 0,fcf = 0,kl = 0), {k\ = M u k\ = 0,Jfcf = 0,Jfe| = 0), (Jfe? = M 1; A;| = M 2 ,fc§ = 
0, fe| = 0), . . .. The coefficients multiplying the poles of Eq. (I2.39P are obtained with 30000 
digits of precision at least (at the FORTRAN level for non-singular points, the proof works at 
the working precision level). The transversality property for the on-shell massless particles 
must be applied also in this case. For the penline and boxline contributions, similar relations 
are found. In addition, for them, we have checked the factorization formula analytically for 
all possible kinematic configurations, once the on-shellness of the gluon, p 2 = 0, is imposed. 

As in the pure abelian case, the factorization of the divergences, for the sum of the Ai^ 
terms, against the born amplitude at the numerical and analytical level already provides a 
strong check of the correctness of the result since Ai and Ai^ have the same analytical 
structure, Eq. (12. lip . The cancellation of the renormalization scale variable, /x, in Aiv, 
Eq. (I2.4ip . provides an additional strong check (see Table [T71 for hexagons). Additionally, 
similar to the pure abelian case, we have implemented Ward identity tests for the virtual 
corrections in FORTRAN and Mathematica at different levels of complexity: 

1) At the level of single diagrams. 

2) For the different topologies M 9 Vl T , M 9 ViV2 T and M VlV2 y 3T . 

3) At the level of gauge invariant amplitudes for one gluon emission. 

4) Subset of amplitudes invariant for a specific replacement (e{pk) — e k Pk)- 

5) Specific contractions that make the contributions to vanish. 

1) and 2) under the replacement of the polarization vector by its momentum (e(pk) = — > 
Pk) for the EW vector bosons, relations among the different diagrams/contributions can be 
obtained. Considering the hexagon of Figure dU 

fjC^^p^p^p^p^p^ps) = ( 2 - 45 ) 

d'q 1 1 1 1 1 1 rra I _ \ 

(2tt) q! q 2 (q - p 6 ) 2 $ + f 1A $ + f u i + f 12 $ + ti 
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with V"(q — Pq,Pq, q), the triple gluon vertex. Note that we have kept the complete depen- 
dency of the momenta for clarity despite the fact that we can use momentum conservation 
to eliminate one of them. Then, contracting one of the open indices by the corresponding 
momentum and expressing the contracted gamma matrix as the difference of two adjacent 
fermionic propagators, we find the following relations, 

(Pl,P2,P3,P4,P5,P6) = ^ 3/14At6 (Pl,P2+P3,P4,P5,P6) ~ ^^^(Pl + P2,P3,P4,P5,Pg) 

P3 3 - 7r ^3M4w(Pl'^2,P3,P4,P5,P6) = £f t * 2 ^ 6 (Pl,P2,P3 + P4,P5,P&) ~ £f 2 ^ 6 (Pl,P2 + P3,PA,P5,Pg) 

P4 4:F C^^ (i (Pl^P2,P3,P4,P5,P6) = £/£ ww (PliP2>P3iP4+P5 ) l»6) ~ £f 2 ^ 3 ^(pi,P2,P3 + P4,P5,Pg), 

(2.46) 

where S^ 1I12I13 represents the pentagon diagram of Figure HU defined similarly as Eq. ( 12. 45ft . 
For the contributions, Vi G (W ± ,Z, 7) is considered such that only one global color factor 
appears for all diagrams of Figures (TOj [Tl] and [12j Then, factorizing out the polarization 
vectors and couplings, 

./<...V„,>1, • • • iPn+2) = gX Xf ■ ■ ■ 9r nf 9oet?M ■ ■ ■ 4:(Pn)e^ +2 (p n+2 )^^... MnM „ +2)r (Pl, • • • ,Pn +2 ) 

(2.47) 

the same relation holds for the hexlineNoAbe and penlineNoAbe contributions under the 
replacements, 

J9* M 9 * 

Zfamm ~^ M 9 ^.^^. (2.48) 

Similar relations are obtained between the penlineNoAbe and boxlineNoAbe contributions. 
Finally, the EW Ward Identity for the boxlineNoAbe routine returns directly zero. 

3) To check gauge invariance for EW boson production with an additional jet and due to 
the fact that QCD mixes, for example, the topologies of the Figures [6] and H2J we have to 
consider the abelian type contributions of the previous section, including the cross terms, 
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Figure 15: Direct and cross terms of Eq. ( 12 . 2 1 [) and Eq. ( 12 . 3 7|) forming a gauge invariant 
quantity. 

along with the non-abelian type contributions (see Figure [15] for the convention of momenta 
and to illustration). For the different topologies, it is satisfied that 

Pf (M/ 12fl3 , T (p 1 ,p 2 ,P3) + M fl3fl2 , T (pi,p 3 ,p 2 ) + M^^r (pi, P2, PA, Pz)) = 

pf (H2Ww,r(Pi,P2,P3,P4) + 5 other + Mf l2 ^ fl4 , T (p 1 ,p 2 ,P3,P5, Pi) + 1 other) = 

pT (>W3«m 5 ,t(Pi,P2,P3,P4,P5) + 23 other + Mf l3fl2ft4fi5 , T (pi,P3,P2,P4„P6,P5) + 5 other) = 0. 

(2.49) 

The abelian contributions with an emitted gluon give rise to two different color structures, 
Fig. ED The non-abelian type contributions are proportional to Ca and their contribution 
cancel in the gauge test Eq. (I2.49j) against the abelian type ones proportional to Ca- The Cp 
terms from the abelian contribution cancel against themselves as discussed in the previous 
section, therefore, each one of the equalities of Eq. ( I2.49P allow us to check on the one hand 
the abelian type contributions and, on the other hand, the interplay between the abelian 
and non-abelian contributions. 

4) For penlineNoAbe and hexlineNoAbe contributions, there are as well "gauge invariant" 
subsets for an specific replacement. For example, for the hexline and hexlineNoAbe contribu- 
tions with gluon momentum, p 5 , an invariant subset under the replacement e(ps) = — > p^ 
is obtained by permuting the position of the momentum and the corresponding Lorentz 

6 Note that for processes with W vector bosons not all the permutations are physically allowed. 
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index and keeping the relative order of the other vector bosons fixed, i.e., 

(M„ WIM ^ iT (p U p-2, Pa, Pi, Ps) + M^ WWIMjT (Pl,P2, Pa, P5, Pi) + 

+M^ 5IX3lMijT (p 1 ,p 2 , Ps, Pa, Pa) + -M |46tfaWM4 , T (pi, p B , P2, Pa, Pa) + 

+- M M2M3M4M5,r(Pl>P2,P3,P4,P6,P5)) =0. (2.50) 

Note also that the divergences of this combination also factorize to the born amplitude, 

Eq. 

5) In addition, we have checked analytically and later on numerically that there are some 
combinations of replacements for the polarization vectors which make the contributions to 
vanish. 

bf) {M^ T ( Pl ,P2,P3,P4)) = 0, 

(pT(P2+P 3 r 3 \\(P2+Par 2 pT) (M^ 3fM5tT (pi,P2,Pa,P4,p 5 )) = 0, (2.51) 

P' 2 L2 {P2+Pa+PAY 3 PT (Mf 12 ^ IM4 ^ T ( PllP 2i P 3,P4,P5,P6)) = 0. 

Note that the first line represents the EW Ward identity for the boxlineNoAbe, equivalent 
to Eq. (12.461) . which is zero as mentioned previously. 

3 Numerical Instabilities and Timing 

To study the stability and timing of the contributions computed here, we have used 5 • 10 5 
cut-accepted points for the vector boson fusion process qq —> qqW + W~ at LO generated 
with VBFNLO [89j applying the cuts, 

p Tj > 20 Gev, \ Vj \ < 4.5, A Vjj > 4, (3.1) 

PTj and \yj\ are the transverse momentum and rapidity of the jets, respectively. Ayjj is 
the tagging jet rapidity separation. Furthermore, we require the two tagging jets to lay in 
opposite detector hemispheres, 

Vji x y j2 < 0, (3.2) 

with an invariant mass of Mjj > 600GeV. For the hexline(NoAbe) contribution, we assign 
the momenta of the two outgoing jets to two of the vector bosons emitted from the quark 
line in Fig. [6] (Fig. [T2l . The corresponding polarization vectors are constructed following 
the conventions of Ref. [92] (Eq. (A. 11) of that paper). For the boxline(NoAbe) and pen- 
line (No Abe) contributions, the momenta of the external particles are combined to obtain 
2—7-3 and 2—7-2 kinematics. Different combinations and/or permutations yield similar 
results, we show here the more unstable observed. 

The contributions are written automatically into FORTRAN modular routines in which 
some flags are incorporated to take advantage of the structure of the result. Note that the 
information of the spinor helicity is contained in the matrix elements, SM i r (Eq. (12.81) ). 
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% CPU/Total 


« 1% 


« 30% 


« 70% 



Table 1: CPU time to revaluate different parts of the code under changes of spinor helicities 
and polarization vectors for the hexline contribution. 

which are computed following the helicity method. Thus, calling the routines with a differ- 
ent spinor helicity only requires to recalculate the reduced set of standard matrix elements 
SM ir , which represents less than 1 % of the total CPU time for the hexline (No Abe) con- 
tributions. Additionally, for changes on the vector bosons due to different helicities or to 
the implementation of gauge tests (e(p n ) —> p n ), only the Flj and SMj )T of Eq. (I2.8P should 
be recalculated, but not the Fj, from Eq. (12.91) . In Tabled! as an example, the summary of 
the CPU time spent to revaluate the different parts of the code, as well as the part of the 
code that has to be revaluated under changes of spinor helicities and polarization vectors, 
are shown for the hexline contribution. 

The organization of the code in dependent /independent loop integral parts is of advantage 
since additional callings of the routines for making Ward identity gauge tests are obtained 
at a lower CPU time cost. To control the stability of our routines for the finite terms of each 
of the contributions (boxline(NoAbe), penline(NoAbe), hexline (No Abe)), for every phase 
space point, all the possible Ward identity gauge tests (e(pj) — > Pj) of Eqs. (12.26112. 46p are 
performed. Except for the boxlineNoAbe, which vanishes for the EW Ward identity gauge 
test of Eq. (I2.46p . and Eq. ( I2.49[) is applied with the replacement e(p g ) —> p g . These equations 
can be rewritten in the more convenient form "abs(a/b)-l=0" such that we obtain normalized 
results which allow us to test the accuracy of our numerical zeros, subject that there is not 
any numerical cancellation in the determination of the "a" and "b" quantity. In that case, 
we use "a-b=0" . In practice, we take the worse normalized value of the Ward identity gauge 
test and compare it with a minimum accepted accuracy, i.e., abs(a/b)-l < "accuracy". We 
define as unstable point the one that does not satisfy the Ward identity gauge test at a 
given accuracy. The factorization proof of the divergences to control the stability is not 
used since they need different input integrals (Eq. (l2.5p ) and the complete routine has to 
be revaluated. The use of the Ward identities gauge test for the finite contributions has 
the additional advantage that the same integrals and tensor coefficients, stored in the Fj 
functions (Eq. (l2.9p ). needed to provide the finite result, Ai v , are used. Moreover, since the 
rational terms and the poles will factorize against the born amplitude, we will not revaluate 
the routines for the poles (Eqs. (12.11112. 12ft ). but, instead, in practical implementations, we 
will use the analytical form which only requires to evaluate the born amplitude. 
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The fraction of unstable points for the different contributions depending on the accuracy of 
the Ward identity gauge tests can be seen in Tables [2] and [31 For the determination of 
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Table 2: Fraction of unstable points, out of the sample of 5 • 10 5 events, depending on the 
accuracy of the gauge test for the abelian contributions Eq. (12.171) . Eq. (I2.19P and Eq. (I2.2ip . 
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Table 3: Fraction of unstable points, out of the sample of 5 ■ 10 5 events, depending on 
the accuracy of the gauge test for the non-abelian contributions Eq. (12.331) . Eq. (I2.35P and 
Eq. (E37|). 

the tensor coefficients, we have applied the Passarino-Veltman tensor decomposition method 
following Eqs. ( 1A.7[) and ( 1A.8I) up to the box level, and the Denner-Dittmaier method for 



pentagons and hexagons applying Eqs. ( 1A.13I) and ( 1A.14j) . To deal with these unstable 
points, some knowledge about the origin of these instabilities is required. It is known that 
small Gram determinants appearing in C and D functions and small Cayley determinants 
for E and F functions result in a loss of precision in the determination of the tensor coeffi- 
cient integrals, therefore, in the amplitudes. To solve this problem within double precision 
accuracy, we try to improve the determination of the tensor coefficients up to the box level. 
First, using a fast implementation of the LU decomposition method which avoids the explicit 
calculation of inverse Gram determinants of Eq. (I A. 71) by solving numerically a system of 
linear equations, Eq. flA.llj) . It turns out that this procedure is more stable close to sin- 
gular regions and reduces considerably the number of identified instabilities without CPU 
penalty, Tabs. H] and El In addition, we apply also special tensor reduction routines for small 
Gram determinants in analogy to Ref. |63j using Eq. ( 1A.12I) . These routines are switched on 
whenever a cancellation in the Gram determinant defined by, 

a ^ Xln,,,i n =l e h...i n Pl ■ Ph ■ ■ - Pn ■ Pi n , . 

^ — r, s-^n i j") K^.O) 

^ L~ii\...i n = \ \^i\—inV\ ' Ph - ' 'Pn' Pi„ \ 

is larger than a given cut off. The milder the cancellation the larger the number of terms in 
the expansion that has to be included to provide a good accuracy. For (Dy) functions, 
we include tensor integrals up to rank 9 (7) which guarantees approximately an error of 
0(A 7 ) (0(A 4 )) for the (Aj) tensor coefficients of rank 2 (rank 3) 0. Such that, e.g. the 



7 For D functions, cancellations in more that one sub-determinant can appear. Nevertheless, the estimate 
seems to work. 
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Table 4: Fraction of unstable points, out of the sample of 5 • 10 5 events, depending on the 
accuracy of the gauge test for the abelian contributions. 
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Table 5: Fraction of unstable points, out of the sample of 5 • 10 5 events, depending on the 
accuracy of the gauge test for the non-abelian contributions. 

cancellation of two digits in the Gram determinant typically results in a precision of 14 (8) 
digits for the (-D^) tensor coefficients of rank 2 (rank 3). Cancellations of two to three 
digits, A = 10~ 2 ' -3 , in the Gram determinant of the D functions represent the borderline, in 
our present set up, for which the special tensor reduction routines for Dij functions provide 
better results than the ones obtained with the LU decomposition. Practically, below the cut 
off, we compare the precision provided by the two methods using sum rules and apply the 
most precise one. For up to penline(NoAbe), the combined method works well leaving the 
remaining instabilities below the per mill level for an accuracy of the gauge test of 10~ 4 . 
Note, however, that for the boxlineNoAbe and the penlineNoAbe routines the behavior is 
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Table 6: Fraction of unstable points, out of the sample of 5 • 10 5 events, depending on the 
accuracy of the gauge test for Dble and QUAD precision for the abelian contributions. 

worse than for the abelian contributions. In these cases, cancellations in the calculation of 
Cayley determinants in C and/or D functions also take place and other special routines have 
to be applied, Ref. [63J. 

For the hexagons, the remaining instabilities are still considerable and, generically, cancel- 
lations in Cayley determinants of the E and F functions are also present. At this stage, the 
simplest and fastest way of rescuing the remaining instabilities is calling the contributions 
with quadruple (QUAD) precision. However, quadruple precision is 20 times slower than 
double (Dble) precision, thus, revaluating the identified unstable points, 5% of the points 
for the hexline contribution (the most time consuming and most unstable object) for an 
accuracy of the Ward identity gauge test of 10~ 6 , results in an addition of 100% CPU time. 
A better approach which reduces the slowing factor of QUAD precision consists in applying 
QUAD precision to compute the input scalar and the tensor coefficient integrals. The tensor 
decomposition routines and scalar integrals amount a fraction of the total CPU time of the 
hexline (No Abe) routines and evaluating these in QUAD precision results in routines only 
a factor 4 slower than the original one. Therefore, for instabilities at the 5% level, only 
an additional 20% CPU time is added, in contrast to the 100% when applying quadruple 
precision to the whole routine. After this procedure, the instabilities are reduced for an 
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accuracy of 10~ 6 approximately a 300 factor for the hexagon contributions, leaving the in- 
stabilities at the ~ 0.15%o level, confirming that the lost of accuracy is at the level of the 
tensor integral coefficient determination due to cancellations in the Gram and/or Cayley 
determinants. For this reduced set of points, ~ 0.15%o, Tab. [HI we can call the routines with 
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Table 7: Fraction of unstable points, out of the sample of 5 ■ 10 5 events, depending on the 
accuracy of the gauge test for Dble and QUAD precision for the non-abelian contributions. 

full QUAD precision with only an additional 5 %o CPU time, reducing the instabilities by 
an order of magnitude. In Tables M arid [7J the fraction of unstable points after these steps 
depending on the accuracy demanded for the Ward identity gauge test are shown for the 
abelian and the non-abelian contributions, respectively, for two setups, with and without 
previously applying the special routines for small Gram determinants. The non-existence of 
instabilities for an accuracy of the Ward identity gauge test of 10 -4 , both for the hexline and 
for the hexlineNoAbe routines, without using the special routines for small Gram determi- 
nants in double precision, suggests that an increase of precision would be enough to rescue 
these points, and the Landau singularities are not present for cut accepted points or they are 
the sub-million level for the kinematics of EW pp — » W + W~jj + X production and given 
the cuts of Eqs. (I3.1|3.2p . Once, we are forced to use QUAD precision for the hexagons, 
the use of the special routines does not bring too much of an improvement concerning the 
final number of instabilities but it introduces a delaying factor. To reduce this factor, we 
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only switch on the special routines whenever the cancellation of the Gram determinant are 
severe, larger than 3 digits. This procedure is used to compute the timing of the routines in 
Tables E] and In these Tables, one finds the average CPU time for the evaluation of the 
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Table 8: Average CPU time per point for the abelian contributions including additional time 
spent performing gauge tests and rescuing unstable points. 

abelian and non-abelian type contributions, respectively, including a detailed description of 
how this time is distributed for the calculation of the Ward identity gauge tests, different 
helicities and the revaluation of unstable points with QUAD precision using a single core of 
an Intel(R) Core(TM)2 Quad CPU Q6600 @ 2.40GHz computer with the Intel FORTRAN 
compiler. We observe that most of the time of the hexline (No Abe) routine is spent in ap- 
plying the Ward identity gauge tests despite the fact that not all of the routine should be 
revaluated, Tab. [TJ We can reduce this slowing factor without affecting statistically the iden- 
tified number of instabilities given in Tables |6] and [7] by applying instead of the four (three) 
gauge tests for the hexline (No Abe) contributions only one selected randomly. The total CPU 
time in parenthesis in Tables M and M is computed applying this procedure. Moreover, the 
time shown is referred to the calculation of the finite pieces Ai v since we will assume that 
the divergent pieces and_A/^ forcpmplete processes are known analytically, such that, we do 
not have to revaluate Ai v and M v for the divergent part of the input functions. 

In Fig. [16J the relative accuracy of the double precision result, eo, defined as the absolute 
value of the difference between the double and quadruple precision divided by the quadruple 
precision result is plotted. We also plot represented by the "QUAD's" label, the difference 
between the amplitudes evaluated with full QUAD precision and with QUAD only applied to 
the basis of scalar integrals and tensor integral routines dotted-dashed line. The effect of the 
rescue system, i.e., the double result is set to the QUAD precision when the Ward identity 
gauge test for the amplitudes fails at the accuracy 10 - ^ are described by the "Impx" lines. 
In the right panels, one can see the critical region, where double precision is not accurate, 
and the effect of the rescue system for different values of the Ward identity gauge test. One 
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Table 9: Average CPU time per point for the non-abelian contributions including additional 
time spent performing gauge tests and rescuing unstable points. 

observes that the choice 10~ 4 (Imp 4 ) gives 3 digits of precision at the per mill level with 
non-identified left instabilities at this level of accuracy, Tables |6] and [71 From the plot, we 
can conclude that the Ward identity gauge test together with the use of QUAD precision is 
an efficient system to control the accuracy of our results. 



4 Conclusions 

In this paper, some of the NLO QCD one-loop amplitudes contributing to the process pp — > 
Wjj + X have been computed. It has been shown that the numerical instabilities due to 
the presence of small Gram and Cayley determinants are under control. The use of the LU 
decomposition method for C and D functions reduces considerably the fraction of unstable 
points at the double precision level (without CPU penalty). This improvement together with 
the use of special routines for small Gram determinants is enough to reduce the instabilities 
for pentagons well bellow the per mill level for an accuracy of the Ward Identity gauge test of 
10~ 4 . Using quadruple precision in two steps, the instabilities for the hexagons for a precision 
of 10~ 6 in the Ward identity gauge tests are reduced up to the 0.03%o level. To decrease 
further this number the use of higher precision can be safely applied since for an accuracy 
of 10~ 4 the identified instabilities disappear even without previously using special routines 
for small Gram determinants. This suggests that the presence of Landau singularities which 
translate into exactly vanishing Gram determinants, for accepted cut points, is small for the 
kinematics of EW pp — > WWjj + X production, given the cuts of Eqs. (13. ip and (13. 2p . 
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Figure 16: The normalized difference between Double precision and Quadruple precision 
for the hexline (top) and the hexlineNoAbe (bottom) contributions is represented by the 
"Raw" line. "QUAD's" represents the difference between the amplitudes evaluated with 
QUAD precision and with QUAD only applied to the basis of scalar and tensor integral 
routines. "Impx" represents the effect of the rescue system activated for an accuracy of the 
Ward identity gauge test of 10~ x . The right panels describe the critical region where double 
precision is not accurate and the effect of the rescue system for different values of the Ward 
identity gauge test. 



It has been shown factorization proofs of the infrared divergences of the bosonic contributions 
of WW and VWj production with V G (W, Z, 7). We have also given the master 
equations for the evaluation of the tensor coefficient integrals in the external convention 
including those needed for small Gram determinants. 

The CPU time of the contributions involving hexagons, which rounds the milliseconds, are 



33 



competitive with other more sophisticated methods. We plan to use these tools not only to 
compute diboson plus two jet production at hadron colliders at NLO QCD, but also other 
interesting 2 — > 4 processes like Wj'j + j [75J. 
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A Tensor Decomposition Master Equations 

For the derivation of the Master Equations presented in this Appendix, we have followed 
closely Ref. [63J. The difference resides in the convention used to derive them. Meanwhile 
Ref. [B3] uses the internal propagator momenta, q k , we use the external momenta, p k - The 
conversion from one notation to the other is trivial. Nevertheless, we have derived the 
recursion relations in this notation, first, to avoid loss of accuracy in the conversion for 
critical points, second, the conversion for pentagons and hexagons of high rank starts to be 
lengthy and time consuming. 

The n-point one loop tensor/scalar integrals can be written as 

l£"**(p 1 ,...,p N - 1 ,m ,...,m N - 1 ) = c r (ji 2 ) / d D l ^ - - , P>0 (A.l) 

with 

^ » = km^m • (A - 2) 

and 

N k = {l + q k f-ml + it= i^ + Y^P^ -ml+ie, k = 0, . . . , N-l, q = p = 0. (A.3) 

Although the introduction of external momenta may seem an additional complication, if 
one carelessly uses the internal propagator notation but uses external momenta in their 
expressions, the number of terms grow factorially. As an example, we present the tensor 
integral of rank 1 for a box integral. With the external momenta convention, the tensor 
integral is written in terms of their tensor integral coefficients as, 

If (pi,Pa,p 3 ) = ^iPf + D^i 1 + Ajpf- (A.4) 
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If we use the internal momenta notation, IM, for the same integral and translate it to the 
external momenta notation to apply equation of motions or transversality properties, the 
number of terms increases from 3 to 6: 

It 12, fc) = D[ M qi + Di M q 2 + Dj% = D™tf + D 2 /M ( Pl + D{ M { Pl 

(A.5) 

For Hexagons of rank 5, this means passing from 4500 to 200000 terms. Although, this 
problem can be solved trivially, e.g., applying conservation of momentum, some care has to 
be taken to avoid exceeding the memory capacities of the computer at intermediate stages. 

Generally, the tensor integrals are written in terms of its tensor integral coefficients using 
Lorentz-covariant structures, 

[P/2] N-l 

= £ £ {^^•••rfw^f...0 > „ + ,., P (A-6) 

n=0 t2t»+i,...,tp=l „ 



2n 



where the curly braces are defined similarly as in Ref. |63j . 

Following Ref. [63], we find for the external momenta convention up to the pentagon level, 
the following Master equations for the tensor coefficient integrals: 

N-l / P 

jN _ V^/7(w-m(-i) [ o(N-i) _ 2 Vr) ■ T N - 

1 i 1 ...i P - / Jim I °ni 2 ...ip ~ Z_^, '■'>■'■ ... 

n=l \ r=2 

N < 5, 1 < i n < N - 1 (A.7) 



ll 3 .., P = 2{D + l p _ N) ((-l)tW--^).)g(0) + 2mlll_ ip+ (A.8) 



N-l 



+ £ ( r « " r n-i)I^3.., P ), N<5,i 3 ...i P = 0,...N-l 



71=1 



(Z^ 1 ^) < - ^ is the inverse Gram determinant which is built with external momenta, Z^ 
2(pi ■ pj) and r n = (j2" =1 pj) - m*. Sf+] ip is defined by 



dij = 1 - Sij, is given by, 



A - , (a.9) 



and j in I^ij) means that the Nj propagator has been canceled, that is, the tensor co- 



efficients for the tensor integral with the Nj propagator removed from For the LU 
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implementation, the inversion of the Matrix in Eq. (IA.7j) is not done analytically, instead, a 
system of equations for the tensor integrals given by, 

p 

(Z (N 1) )i in lZ..i P = S ni 2 ...i p -' 2 ^2 5 ni r lQ 0i2 _l r ip , (A- 11 ) 

r=2 

is solved using a fast implementation of partial pivoting of the LU decomposition method. 
For small Gram determinants, we define in analogy to Ref. [63], 

s^hp = s N(i2h . ■ .^ p)j ./ ( T 2) ,.. Mi (i) - (-i)^'^-^-^^'-^^),-,..^^,^ - 1) 

«U = 2(-l)^.o^^»o)/J )o ... (!p) jO) + 2mo /^, ip)o) 

(A.12) 

and the recursion relations are obtained directly from Eqs. (5.41-5.48) of Ref. [63] using 
external momenta in all (sub) determinants and the replacement f n — > (r n — r n _x). We have 
implemented these recursion relations up to rank 9 for the CV,, and up to rank 7 for the -Dis- 
functions 

For the pentagons, we have followed closely section 6 of Ref. [63] and derived their master 
formulas Eqs. (6.12,6.13) for the external momenta convention. Although, the cancellation 
that takes place to get their simple equations are different, the final master integrals are 
quite similar: 

4 

det(X^)E kH ... tp = Y^xtl {Kn...S ip nD {h)n ... iip)n (n)- (A.13) 

n=l 

-(-l) (5o(il) ^ 1+ - +5o£ - ) ^ l) ^(n) n -,..(ip) n - 1 (« - 1)) - A2 ) (-l) (5o(ll)o+ - +5o( - )o) J D(n)o-...(.p)o(0) 

4 P 

~ 2 S S X (kn)(0ir) • • • K-mK+m ■ ■ ■ ^pn jD 00(ii)„...(iOn...(ip)„( n ) _ 



n=l r=l 



V 1 00(il) n _l.. ■(^r)?i — 1 ••• 



fc = l,...,4, P<4, 



4 

rfet(X (4) )E 00 i 2 ...ip = y^nlJ (<^2« • • • ^pn-D 0(i 2 ) n ...(i P )„(™)- 
n=l 

-(-l)^ 2 ^^-^ ^^^^)^,..^)^^ - 1)) , P < 4, (A.14) 

where X( 4 ) and related quantities are defined/obtained from Eq. (2.23) in Ref. [63J making 
the replacement /„ — > (r n — r n _i) and using directly external momenta, — > 

8 The UV part of the integrals are not given explicitly here. They can be obtained by contacting the 
author. 
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For the hexagons, we have followed section 7 of Ref. [63]. The master equations with ex- 
ternal momenta convention are obtained from their Eqs. (7.13,7.16-7.18) by making the 
replacement: 

E i2i3 .., P (0) -»> (-l) ( ^^ 1+ "^ ) «- l) %)n-i...(*,.)»-i(n " 1), ( A -!5) 

using directly external momenta, —±Pk- We have checked the derivation of all the ten- 
sor coefficients up to Hexagons of rank 5 by contracting the tensor integrals with external 
momenta. Thus, we can relate all the tensor coefficients of order N and rank P to tensor 
integrals of order N and N — 1 and rank P — 1 : 

N-1 

r> r> n, TUl-PP — V~" 7 7 7, T N 

Ft^iFj^ • • ■ yk(ip 1 JS[ / J /J ll\ £J ji% ■ ■ ■ ^kip L i x ...i p 

i\...ip 

Pi* ■ -^(CrW - (< - 1) - fa - r.-i)/^-^). (A.16) 

For hexagons of rank 5, up to 166 different combinations (not all of them independent since 
only four momenta out of six are linearly independent) where built and checked at the 
FORTRAN and Mathematica level. The implementation of these routines together with 
the massless input integrals in Mathematica have proved to be quite useful since we can 
overcome the loss of accuracy due to small Gram determinants by increasing the working 
precision, and check accurately the special routines for small Gram determinants, the tensor 
reductions for pentagons a la Passarino Veltman vs a la Denner-Dittmaier as well as the 
complete numerical result of the contributions. 



B Benchmark numbers for the hexagon contributions 



In this section, we provide numbers for the hexline and the hexlineNoAbe contributions. 
The set of momenta used for q(pi) q(pz ) Vi V2 V3 g — > is chosen to be: 



Pi 

P2 
PVi 
PV 2 

Pv 3 
Pa 



2210.591640000411, 0, 0, 2210.591640000411), 
410.6465697388802, 0, 0, -410.6465697388802), 

-1644.598252136518, -186.7167811992445, 14.28499809343437, -1633.902137026130), 
-266.9005208707261, -59.67643702980636, -176.1526035584612, -173.8364669053697), 
-402.5905961231937, 166.3844370751274, 218.6277243718783, -283.2688206929862), 
-307.1488406088535, 80.00878115392345, -56.76011890685155, 291.0623543629552), 

(B.1) 



with the notation p = (p° , p 1 , p 2 , p 3 ) and all the components given in GeV. This represents a 
cut-accepted point for the EW qq —> W + W~qq process with a "mild" behavior, with "mild" 
meaning that there are not small Gram/Cayley determinants appearing up to the hexagon 
level, therefore the double precision routines provide accurate results. The corresponding 
polarization vectors are constructed following the conventions of Ref. [32], Eq. (A. 11). First, 
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we give results for each of the diagrams that constitute the hexline contribution for a par- 
ticular permutation of the electroweak vector bosons, M.ViV 2 v 3 g,r, with r = —1, in Table ITUl 
and Table [UJ using the decomposition, 

Re(M*-M* VlV2V3g , T ) _ gl " (C ^\ 

2 \M?M?*\ " (4^ (47r)r(1 + e)( ^ o) 2-,{—f^) x 

(n) C (l),ViV 2 V3g , (-D ~ 4) ( n ) (-P - 4) d( (i),v 1 v 2 V3g \ 

C (0),V 1 V 2 V :i g + e + £ 2 + _ 2e a (0),ViV 2 V3 9 + _ 2e £ I ' 



correspondingly to Eq. ( I2.2ip . M.^. = A4y lV2 y 3gT , i.e. the born amplitude for the specific 
order of vector bosons. The (47r) e r(l + e)(/io) -2<E factor comes from the definition for the 
scalar and tensor integrals, Eq. (1A.1I) . We can implement some tests to these numbers. 





(n) 

C (0),ViV 2 V 3 g 


(n) 

C (l),ViV 2 V 3 g 


(n) 

C (2),V 1 V 2 V 3 g 


n= 1 


-0.8495224242548815E+02 


-0.2718191945175511E+02 


-0.3999999999999874E+01 


n= 2 


0.73737152265451 16E+02 


0.1475364608414150E+02 


-0.4823858365463354E-13 


n= 3 


0.7720679379271682E+02 


0. 1032635756877249E+02 


0.4778146726503880E-13 


n= 4 


-0.7023321057368833E+01 


-0.2247107183707837E+01 


-0.2411275830791782E-14 


n= 5 


0.5075981735087911E+01 


0.0000000000000000E+00 


0.0000000000000000E+00 


n= 6 


-0.1241513230258542E+03 


-0.1483523289520358E+02 


0.3441661498683095E-13 


n= 7 


0.1268199899143818E+02 


0.2000000000000156E+01 


0.1502284155171117E-26 


n= 8 


0.9572211011177236E+01 


0.2000000000000001E+01 


0.0000000000000000E+00 


n= 9 


0.2215410796479108E+02 


0.2000000000000001E+01 


0.0000000000000000E+00 


n= 10 


0.1060058346584089E+03 


0.1318425587774655E+02 


-0.1013265897019504E-13 


n=ll 


-0.1268199899143716E+02 


-0.2000000000000000E+01 


0.0000000000000000E+00 


n=12 


-0.1024742028732465E+02 


-0.2000000000000000E+01 


0.0000000000000000E+00 


n=13 


-0.1523236618746480E+02 


-0.2000000000000000E+01 


0.0000000000000000E+00 



Table 10: Contributions of the hexline diagrams. c u)vvvo are defined in Eq. (1B.2j) . 



Assuming a general common factor, Cp, we can check the value of the sum of c^j lVl v 2 v s g ^ 
testing the factorized scale energy independence, similarly as Eq. (12.411) through, 

E Ws>°) - 2Re (f(s, t, it, // , 1, 0, 0)*) £ F(/i ) (B.3) 

with f(s, u, t, /i , 1, 0, 0) given by Eq. (I2.42p with s = (pi +P2) 2 ,t = (pi + p g ) 2 and u = 

(p 2 + P g ) 2 - In Table [T2j we show the values of the sum of c^} Vl v 2 v 3 g anc ^ ^( s ' M > Z 10 ' ^ ^> ^) 
for two different sets of factorization energy scales, for /i = s, corresponding to the values 
given in Table [TDJ, and for /i = lGeV. One can see that Eq. (IB. 3D is satisfied with an 
accuracy of 12 digits. Additionally, for this specific phase space point, the Ward identities of 
Eq. (12.261) for the complete contributions are satisfied at the 12 digit level. For the divergent 
contributions and the rational terms, the factorization against the born amplitude is shown 
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Jn) 

a (0),ViV 2 V 3 g,T 


u (l),ViV 2 V 3 g,T 


n= 1 


0.8702382558218599E-12 


-0.8613060482805703E-13 


n= 2 


0.1709040047358825E-13 


-0.8545281855364259E-15 


n= 3 


0.3736943504519546E-13 


-0.4397635645750366E-14 


n= 4 


-0.2396002150441086E-13 


0.1363187596482680E-15 


n= 5 


0.0000000000000000E+00 


0.0000000000000000E+00 


n= 6 


-0.3591236576358064E-12 


-0.4690957329075479E-14 


n= 7 


-0.2000000000000079E+01 


0.0000000000000000E+00 


n= 8 


-0.2000000000000001E+01 


0.0000000000000000E+00 


n= 9 


-0.2000000000000001E+01 


0.0000000000000000E+00 


n= 10 


-0.2000000000000000E+01 


-0.6325727406404221E-14 


n= 11 


0.2000000000000001E+01 


0.0000000000000000E+00 


n= 12 


0.2000000000000001E+01 


0.0000000000000000E+00 


n= 13 


0.2000000000000001E+01 


0.0000000000000000E+00 



Table 11: Contributions of the hexline diagrams. d^v^Vsg are defined in Eq. flB.2j) 







Ho = lGeV 


2-u C (0),ViV 2 V 3 9 


0.5214540844413348E+02 


-0.3135489678574997E+03 


2Re(f(s,t,u,ft , 1,0,0)*) 


0.1973920880217872E+02 


-0.3459551674994638E+03 


Eq. (jB.3) 


0.3240619964195476E+02 


0.3240619964196406E+02 


ratio- 1 


-0.2868816295631405E-12 



Table 12: Check of factorization scale energy independence through Eq. flB.3j) for two dif- 
ferent sets of factorization energy scales. 





(n) 

C {l),ViV 2 V 3 g 


H 

C (2),ViV2V39 


d {n) 

U (0),ViV 2 V 3 g 


E 


-0.6000000000005826E+01 


-0.3999999999999853E+01 


-0.1999999999999537E+01 


Exact 


-0.6000000000000000E+01 


-0.4000000000000000E+01 


-0.2000000000000000E+01 


ratio- 1 


0.9710010573371619E-12 


-0.3674838211509268E-13 


-0.4524345419647128E-12 



Table 13: Factorization of the divergent contributions through Eq. flB.4j) against the born 
amplitude and accuracy. 



in Table [13] through, 

E(«) _ a 

C (l),ViV 2 V 3 g - 



(n) 

(0),ViV 2 V 3 g 



(B-4) 

correspondingly to Eq. (I2.44p . Note that the numerical values for dmviViVsg * n Table [TT] are 
below lO^ 13 ). 

In the following, we give results for the hexlineNoAbe contribution for r = — 1. We follow 
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9i{ n ) 

c (o)yiV 2 v 3 


C {l),ViV 2 V 3 


r, ( m\ 

9A n ) 
C {2),V 1 V 2 V 3 


n= 1 


-0.3970812374704999E+02 


-0.1265721313949629E+02 


-0.1813569591709497E+01 


n= 2 


-0.1350281050903405E+01 


-0.7519564890398996E+00 


-0.2022904490007958E+00 


n= 3 


0.2728377087017317E+02 


0.2055535050768235E+01 


-0.5549746075488889E+00 


n= 4 


-0.3264136921557698E+00 


-0.1238348703073990E+00 


-0.1416214405745540E-01 


n= 5 


0.1723360559577099E+01 


0.3328092079219656E+00 


0.0000000000000000E+00 


n= 6 


-0.5023976788785847E+02 


-0.1779284983273073E+02 


-0.2954984559380555E+01 


n= 7 


0.1253736490780382E-02 


-0.2491704328619187E-03 


-0.7861847392712649E-03 


n= 8 


0.4621898462078639E+00 


0.1193463736617884E+00 


0.0000000000000000E+00 


n= 9 


0.1418038518281030E+01 


0.4012945693937532E+00 


0.0000000000000000E+00 


n= 10 


0.1937597641023809E+02 


0.8653232589507967E+01 


0.1540767536436781E+01 



Table 14: Contributions of the hexlineNoAbe diagrams, cfjfv^v-i are defined in Eq. ( IB. 51) . 



the decomposition of Eq. flB.5j) assuming EW production with an additional jet, thus, the 



color factor is Ca for all the diagrams, 

2^=^ = A-^rni + e)M~^C A x (B.5) 



\M?Mf*\ ~ (4tt) 



r 9,(n) o,(n) A \d 9 ' {n) 

c {o)y 1 v 2 v 3 + e + e 2 + _2 e (o)yiV 2 y 3 + _2 e e 
correspondingly to Eq. (I2.37p . with, 

that is, the sum of the four Born amplitudes for a given order of the EW vector bosons. The 
individual contributions are given in Table HH and Table [15j In addition, in Table [16] we 
show the hexline results for the four contributions that are mixed under the QCD group for 
a given permutation of EW vector bosons since they will help us to construct useful checks 
for the hexlineNoAbe contributions. We use the notation of Eq. flB.2j) with M.^ given by 



Eq. (IB.6j) and provide separate numbers for the sum of the two different color structures 



that appear in this case. We can check the finite pieces, cfe , by testing the factorized energy 
scale, /io, independence for the two color factors, Cf and C A through, 



C h 



C (&)yiV 2 V 3 g(^o) + c (0)y 1 V 2 gV 3 (f i o) + c (Q)y 1 gV 2 V 3 g(f i o) + c (0), gVxVaVzgi^)) ~ 

-2Re (/(a, t, u, /i , 1, 0, 0)*) £ F(/i ) (B.7) 



-1 



C A 

~2 ( c (o)yiV 2 V39 + c (o),V!V 2 gV 3 + c (o),v ig v 2 v 3 + c {Q),gViV 2 v 3 ) + X/ 

-2Re (f(s, t, u, /i , 0, 1, 0)*) £ F(/i ) (B.8) 



^(0),ViV 2 V 3 
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j9,( n ) 

a (0),ViV 2 V 3 


U (l),ViV 2 V 3 


n= 1 


-0.14364839545681 11E-11 


0.3386511581831592E-13 


n= 2 


-0.6922590837003719E-18 


0.193912241 1552382E-16 


n= 3 


0.5899544848540121E-13 


0.1263944582753722E-13 


n= 4 


-0.5274108820645013E-16 


-0.1071274072706729E-17 


n= 5 


0.3633510145273396E-19 


0.0000000000000000E+00 


n= 6 


0.5975463630574398E-13 


0.1078373100666381E-14 


n= 7 


-0.7861847392712663E-03 


0.0000000000000000E+00 


n= 8 


-0.1416214405745516E-01 


0.0000000000000000E+00 


n= 9 


-0.1371873399923628E+00 


0.0000000000000000E+00 


n= 10 


-0.8478643312109194E+00 


0.1159003560547673E-14 



Table 15: Contributions of the hexlineNoAbe diagrams, dfy ]y x y 2 y 3g are defined in Eq. ( IB. 51) . 







C F - l/2C A 


c (0),gViV2V3 


0.1186972860337361E+00 


-0.2513048326546336E+00 


c (0),VigV 2 V 3 


-0.3260632362384948E+00 


0.1175575424874967E+01 


c (0),ViV 2 gV 3 


-0.3587963452897030E+01 


0.6010092883764584E+01 


c (0),ViV 2 V 3 g 


0.6748418542481737E+02 


-0.2188150305629134E+02 


c (l),9ViV 2 V 3 


0.3435253401752263E-01 


-0.2242562120604109E-01 


C(l),V ig V 2 V 3 


-0.1745714679327207E+00 


0.9440497420994520E-01 


c (l),ViV 2 gV 3 


-0.9768988228536071E+00 


-0.1710884136252724E+00 


C (l),ViV 2 V 3 g 


0.1091241185529769E+02 


-0.1569618503789184E+02 


C(2),gViV 2 V 3 


0.0000000000000000E+00 


-0.3144738956157922E-02 


C(2),V ig V 2 V 3 


0.0000000000000000E+00 


-0.5664857622981814E-01 


C(2),y 1 y 2g y 3 


0.6661338147750939E-15 


-0.5487493599694454E+00 


C(2),Viy 2 V3g 


-0.4263256414560601E-13 


-0.3391457324843469E+01 


d(0),gViV 2 V 3 


-0.2765908162227770E-14 


0.7861847503483757E-03 


d(0),VigV 2 V 3 


0.4302114220422482E-15 


0.1416214405720360E-01 


d(0),ViV 2 gV 3 


0.4385380947269368E-14 


0.1371873399918132E+00 


d{0),ViV 2 V 3 g 


0.3652633751016765E-13 


0.8478643312106809E+00 



Table 16: Contributions of the hexline diagrams for the different colors and for the four 
permutations for a fixed order of the EW vector bosons. 



similarly to Eq. ( I2.4ip . The prime quantities are built by adding the two numbers of each 
line in Table [TBI For Eq. ( 1B.8j) . C(o) is obtained by taking the number of the second column of 
each line in Table [16j We give the results for two different sets of factorization scale energies 
in Table HTl and the accuracy of it. We obtain 11 digits of precision for the specific phase 
space point. Finally, proceeding in a similar manner from Eqs. (I2.43|2.44p . the factorization 
of the divergence contributions can be checked for the two color factors, the results are 
shown in Table [T81 The value for all of the dm quantities is always below 10~ 13 (not shown). 
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yWo = 1 




Cp color factor 


First line of in Eq. (|B.7|) 


0.4874171644140917E+02 


-0.3169526598603159E+03 


2Re( f(s,t,u,uo, 1,0,0)*) 


0.1973920880217872E+02 


-0.3459551674994638E+03 


Eq. (IB.7) 


0.2900250763923045E+02 


0.2900250763914789E+02 


ratio- 1 


0.2846611835138901E-11 




Ca color factor 


First line of Eq. (1B.8j) 


-0.4135999643699960E+02 


-0.1991502907112327E+03 


2Re(/( S ,t, M ,/i , 0,1,0)*) 


0.7711028150622136E+02 


0.1260325556690930E+03 


Eq. (|R8) 


0.4322385485937497E+02 


0.4322385485945347E+02 


ratio- 1 


-0.1816102823681831E-11 



Table 17: Check of factorization scale energy independence through Eqs. (1B.7IIB.8I) for two 
different sets of factorization energy scales. 



"Exact" for C(i) is given by the 16 first digits of (See Eqs.f l2.43|B.2|B.5j) ): 



2 . Re[ (-lo g (z£) + lo g (^) + lo g (^ 


))*]■ (B-9) 




C(l) 


C(2) 






Cf color factor 


E 


-5.99999999998432 


-3.99999999999893 


-1.00000000001008 


Exact 


-6.00000000000000 


-4.00000000000000 


-1.00000000000000 


ratio- 1 


-2.612909888455306E-012 


-2.668976151198876E-013 


1.008459982188015E-011 




Ca color factor 


E 


-11.8662386614969 


-2.00000000000024 


< E-13 


"Exact" 


-11.8662386614914 


-2.00000000000000 





ratio- 1 


4.600764214046649E-013 


1.187938636348917E-013 





Table 18: Factorization of the divergent contributions against the born amplitude and ac- 
curacy for the two color factors. 



C Color factors 



The color factors for the two permutations of the boxline for one gluon emission are given 
by: 

_ ( -(l) _ L (2) ~~ (3) ~~ 1 y^ F ~ l l lL/ A) 

r v 9 _ r>9 v — r9 v — r v 9 _ (n 1 

L (2) - L (3) _ U (4) - L (4) - 1 y^- 1 ! 
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The color factors for the penline routines for one gluon emission are given by: 

(gv 1 v 2 ) 

ngV\V'l _ ngV\Vl _ ngV\V<2, _ rpairi i in si \ 
L (l) ~ L (2) ~ L (A) - 1 K^F — 



(v ig v 2 ) 



(1) - -(2) - -(4) 

ngV\Vi _ ngV\V2 _ ngV\Vi _ ngV\Vi _ ngV\V<2, _ rpari /p r>\ 

^(3) ~~ U (5) ~~ U (6) ~~ U (7) ~~ U (8) ~~ 1 



nV\gV2 r>V\gVi nV\gVi nV\gVi Tatri -r Icisi \ 

L (l) _ L (2) _C -(3) _C (5) _J {^F-i-f^A), 

pVigV'2 nV\gVi nV\gVi nV\gVi rpari ( ri q\ 

L (4) - L (6) ~ L (7) - L (8) ~ 1 U F: {^^) 



L (l) - L (3) ~ L (6) ~ L K U F-i-/^A), 

nV\Vig _ pV\V2g _ nV\V<2,g _ nV\Vig _ nV\V<2,g _ rpari /p a\ 

^(2) ~ ^(4) ~ U f5) ~~ U m ~~ ^(8) ~ 1 



the same color factors are obtained for the permutations corresponding to V\ -H- V 2 with ]/} 

e (W±,z, 7 ). 

Finally, the color factors for the one gluon emission case for the hexline routine are given by: 

(gVMVs) 

qVMVs = C gViV 2 v s = qv^v, = T a^ Cp _ 1 / 2 C A ) 

L (3) ~ L f5-6) ~ W8-13) F \^-0) 



(V l9 V 2 V 3 



nV^gViVz /7V1PV2V3 rpai/ri ir,si \ 

C (1 _ 5) = C (g) = 1 {C F - L/ZCa) 



'(1-5) - ^(8) 

/7V15V2V3 _ pV\gV2Vs _ rpari iri fi \ 

L (6-7) - C (9-13) F I^-OJ 



/7V1V23V3 _ pV\V2gVz _ pV\V2gVz _ rpalri 1 Icyri \ 

L (l-3) - L (5-6) ~ L (9) _ J ~ V^Aj 

L (4) ~ L (7-8) _C -(10-13) _ 1 °^ 

(W 2 ^3<?) 

/7V1V2V33 _ ^ViV^Vsg _ pV^VzVzg _ pV\V2Vzg _ rpai ri 1 icyri \ 
L (l) _L (3) ~ L (6) ~ L (10) - 1 K^F-l-I^A) 

nV\V2Vzg _ nV\V2Vzg _ nV\V2Vzg _ nV\V2Vzg _ rpari /p o\ 

U (2) — (4-5) — (7-9) — U (ll-13) — 1 \^- ) 

the same color factors are obtained for the six other permutations of the vector bosons Vi, V 2 
and V 3 with Vi E (W ± ,Z, 7 ). 
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